#include "float.hpp"
#include "getresult.hpp"
#include <random>
#include <vector>
typedef union {
    int i;
    float f;
} FL;
typedef union {
    long int i;
    double d;
} DL;

//穷举[start,end]区间内所有的half点，墙砖为double与mpfr比较误差，输出导致最大误差的x和ulp
float HalfFunction::computeHalfToDoubleAllULP(const half &start, const half &end) {
    cout << "正在穷举[" << start << ", " << end << "]所有half点" << endl;
    short hl_tmp;
    double ULP = 0.0, x = 0.0;
    for (short i = *(short *)&start; i <= *(short *)&end; ++i) {
        hl_tmp = i;
        double inputDouble = (double)(*(half *)&hl_tmp);
        double f1 = getHalfToDoubleOfOrigin(*(half *)&hl_tmp);
        double ulp = getULP(inputDouble, f1);
        if (ULP < ulp) {
            ULP = ulp;
            x = inputDouble;
        }
    }
    printf("x = %.16f, floatx = %.8f, maximumULP = %.2lf\n", x, (float)x, ULP);
    cout << endl;
    return (float)x;
}

//计算所有float点强转为double之后的误差
double FloatFunction::computeFloatToDoubleAllULP(const float &start, const float &end) {
    cout << "正在穷举[" << start << ", " << end << "]所有float点" << endl;
    FL fl_start, fl_end, fl_origin, fl_oracle, fl_tmp;
    fl_start.f = start, fl_end.f = end;
    //去除起点从0开始的若规范数，减少时间开销
    if (fl_start.f == 0.0) {
        fl_start.f = 0.01;
    }
    cout << "一共需要穷举" << abs(fl_start.i - fl_end.i) << "个float点" << endl;
    double ULP = 0.0, x = 0.0;
    for (int i = fl_start.i; i <= fl_end.i; ++i) {
        fl_tmp.i = i;
        double inputDouble = (double)fl_tmp.f;
        double f1 = getFloatToDoubleOfOrigin(fl_tmp.f);
        // double ulp = test2(inputDouble, f1);
        double ulp = getULP(inputDouble, f1);
        // printf("%.16lf %.16lf\n", inputDouble, ulp);
        if (ULP < ulp) {
            ULP = ulp;
            x = inputDouble;
        }
    }
    printf("x = %.8f, doublex = %.16lf, maximumULP = %.2lf\n", x, x, ULP);
    cout << endl;
    return x;
}
//计算float点强转为double均匀number个点误差
void FloatFunction::computeFloatToDoubleRandomULP(const float &start, const float &end, const long int &number) {
    cout << "float数据类型[" << start << ", " << end << "]," << number << "均匀采样" << endl; 
    float step = (end - start) / (float)number;
    FL fl_x, fl_origin, fl_oracle;
    double x = 0.0, ULP = 0.0;
    for (int i = 0; i < number; ++i) {
        float input = start + (float)(i * step);
        fl_x.f = input;
        double inputDouble = (double)fl_x.f;
        double f1 = getFloatToDoubleOfOrigin(fl_x.f);
        // double ulp = test2(inputDouble, f1);
        double ulp = getULP(inputDouble, f1);
        if (ULP < ulp) {
            ULP = ulp;
            x = inputDouble;
        }
    }
    printf("x = %.16f, ULP = %.2lf\n", x, ULP);
    cout << endl;
}

void DoubleFunction::computeAllULP(const double &start, const double &end) {
    cout << "正在穷举[" << start << ", " << end << "]所有double点" << endl;
    DL dl_start, dl_end, dl_origin, dl_oracle, dl_tmp;
    dl_start.d = start, dl_end.d = end;
    double ULP = 0.0, x = 0.0;
    for (int i = dl_start.i; i <= dl_end.i; ++i) {
        dl_tmp.i = i;
        double f1 = getDoubleOfOrigin(dl_tmp.d);
        double ulp = getULP(dl_tmp.d, f1);
        if (ULP < ulp) {
            ULP = ulp;
            x = dl_tmp.d;
        }
    }
    printf("x = %.16f, maximumULP = %.2lf\n", x, ULP);
    cout << endl;
}

void DoubleFunction::computeAllRandomULP(const double &start, const double &end, const long int &number) {
    cout << "double数据类型[" << start << ", " << end << "]," << number << "均匀采样" << endl; 
    double step = (end - start) / (double)number;
    DL dl_x, dl_origin, dl_oracle;
    double x = 0.0, ULP = 0.0, origin_relative = 0.0, relative = 0.0;
    for (int i = 0; i < number; ++i) {
        double input = start + (double)(i * step);
        dl_x.d = input;
        double f1 = getDoubleOfOrigin(dl_x.d);
        double ulp = getULP(dl_x.d, f1);
        if (ULP < ulp) {
            ULP = ulp;
            x = dl_x.d;
        }
    }
    origin_relative = getDoubleOfOrigin(x);
    relative = getRelativeError(x, origin_relative);
    printf("x = %.16f, maximumULP = %.2lf, maximumRelative = %e\n", x, ULP, relative);
    cout << endl;
}
//第二种检测方法，操作double的尾数位和增加显著误差判断分支,其中start和end都大于等于0
void DoubleFunction::detectMethod2(const double &start, const double &end) {
    //模拟半精度输入，尾数位保留10位，其余为全置为0
    DL dl_half_start, dl_half_end;
    dl_half_start.d = start;
    // dl_half_start.i = dl_half_start.i & 0B0111111111111111111111000000000000000000000000000000000000000000;
    dl_half_start.i = dl_half_start.i & 0x7FFFFC0000000000;
    dl_half_end.d = end;
    // dl_half_end.i = dl_half_end.i & 0B0111111111111111111111000000000000000000000000000000000000000000;
    dl_half_end.i = dl_half_end.i & 0x7FFFFC0000000000;
    cout << "检测区间[" << dl_half_start.d << ", " << dl_half_end.d << "]" << endl;
    cout << "预处理采样点数：" << (dl_half_end.i - dl_half_start.i) / 0x40000000000 << endl;
    double input_x = 0.0, ULP = 0.0, relative = 0.0;
    for (long int i = dl_half_start.i; i <= dl_half_end.i; i += 0x40000000000) {
        DL dl_input;
        dl_input.i = i;
        double input_one = dl_input.d;
        double origin = getDoubleOfOrigin(input_one);
        double ulp = getULP(input_one, origin);
        if (ULP < ulp) {
            ULP = ulp;
            input_x = input_one;
        }
    }
    double origin_relative = getDoubleOfOrigin(input_x);
    relative = getRelativeError(input_x, origin_relative);
    printf("预处理阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
    //操作double的后10位尾数位得到一个最大误差，变量定义为ULP，对ULP做判断是不是显著误差，设定分支根据ULP选择不同的检测策略
    //显著误差现在暂且手动设置
    if (ULP <= 100) {
        cout << "******表达式无显著误差,采用检测方法1检测误差******" << endl;
        double ULP2 = 0.0, input_x2 = 0.0;
        //改变检测策略，不从half开始，增大尾数位
        DL dl_float_start, dl_float_end;
        dl_float_start.d = start;
        dl_float_start.i = dl_float_start.i & 0x7FFFFFFFE0000000;
        dl_float_end.d = end;
        dl_float_end.i = dl_float_end.i & 0x7FFFFFFFE0000000;
        //采样
        /* long int sample_number = 100000;
        long int step = (dl_float_end.i - dl_float_start.i) / sample_number;
        for (long int i = 0; i <= sample_number; ++i) {
            DL dl_input;
            dl_input.i = dl_float_start.i + i * step;
            double origin = getDoubleOfOrigin(dl_input.d);
            double ulp = getULP(dl_input.d, origin);
            if (ULP2 < ulp) {
                ULP2 = ulp;
                input_x2 = dl_input.d;
            }
        } */
        //模拟退火
        vector<double> vec1 = simulated_annealing(dl_float_start.d, dl_float_end.d);
        input_x2 = vec1[0];
        ULP2 = vec1[1];
        origin_relative = getDoubleOfOrigin(input_x2);
        relative = getRelativeError(input_x2, origin_relative);
        printf("第一阶段x = %.8lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x2, ULP2, relative);
        DL dl_double_start, dl_double_end;
        dl_double_start.d = input_x2;
        dl_double_end.d = input_x2;
        dl_double_start.i = dl_double_start.i & 0x7FFFFFFFE0000000;
        dl_double_end.i = dl_double_end.i | 0x000000001FFFFFFF;

        //采样
        /* step = (dl_double_end.i - dl_double_start.i) / sample_number;
        for (long int i = 0; i <= sample_number; ++i) {
            DL dl_input;
            dl_input.i = dl_double_start.i + i * step;
            double origin = getDoubleOfOrigin(dl_input.d);
            double ulp = getULP(dl_input.d, origin);
            if (ULP2 < ulp) {
                ULP2 = ulp;
                input_x2 = dl_input.d;
            }
        } */
        //模拟退火
        vector<double> vec2 = simulated_annealing(dl_double_start.d, dl_double_end.d);
        input_x2 = vec2[0];
        ULP2 = vec2[1];
        origin_relative = getDoubleOfOrigin(input_x2);
        relative = getRelativeError(input_x2, origin_relative);
        printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x2, ULP2, relative, log2(ULP2 + 1));
    } else {
        //半精度遍历完之后得到一个触发最大误差的输入，对该输入周围的数操作更多的尾数位
        cout << "******表达式存在显著误差,采用检测方法2检测误差******" << endl;
        origin_relative = getDoubleOfOrigin(input_x);
        relative = getRelativeError(input_x, origin_relative);
        printf("第一阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
        double input_second = input_x;
        DL dl_float_start, dl_float_end;
        dl_float_start.d = input_second;
        // dl_float_end.i = dl_float_start.i | 0B0000000000000000000000111111111111100000000000000000000000000000;
        dl_float_end.i = dl_float_start.i | 0x000003FFE0000000;
        cout << "第二阶段采样点数：" << (dl_float_end.i - dl_float_start.i) / (0x20000000) << endl;
        for (long int i = dl_float_start.i; i <= dl_float_end.i; i += 0x20000000) {
            DL dl_input;
            dl_input.i = i;
            double input_two = dl_input.d;
            double origin = getDoubleOfOrigin(input_two);
            double ulp = getULP(input_two, origin);
            if (ULP < ulp) {
                ULP = ulp;
                input_x = input_two;
            }
        }
        origin_relative = getDoubleOfOrigin(input_x);
        relative = getRelativeError(input_x, origin_relative);
        printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
        //此时尾数位增大到float的23位之后，此时获得的出发最大误差的输入为input_x，此时将剩下的全部尾数位不能穷举了，只能采样检测
        //采样数设置为100000个点
        DL dl_double_start, dl_double_end;
        dl_double_start.d = input_x;
        // dl_double_end.i = dl_double_start.i | 0B0000000000000000000000000000000000011111111111111111111111111111;
        dl_double_end.i = dl_double_start.i | 0x000000001FFFFFFF;
        //采样
        /* int sample_number = 100000;
        double step = (dl_double_start.d - dl_double_end.d) / (double)sample_number;
        for (long int i = 0; i <= sample_number; i++) {
            double input_three = dl_double_start.d + i * step;
            double origin = getDoubleOfOrigin(input_three);
            double ulp = getULP(input_three, origin);
            if (ULP < ulp) {
                ULP = ulp;
                input_x = input_three;
            }
        }
        cout << "第三阶段采样点数:100000" << endl; */
        //模拟退火
        vector<double> vec3 = simulated_annealing(dl_double_start.d, dl_double_end.d);
        input_x = vec3[0];
        ULP = vec3[1];
        origin_relative = getDoubleOfOrigin(input_x);
        relative = getRelativeError(input_x, origin_relative);
        printf("第三阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x, ULP, relative, log2(ULP + 1));
    }
}
//start和end都小于等于0的情况
void DoubleFunction::detectMethod2_1(const double &start, const double &end) {
    DL dl_first_start, dl_first_end;
    dl_first_start.d = start, dl_first_end.d = end;
    dl_first_start.i = dl_first_start.i & 0xFFFFFC0000000000;
    dl_first_end.i = dl_first_end.i & 0xFFFFFC0000000000;
    cout << "检测区间[" << dl_first_start.d << ", " << dl_first_end.d << "]" << endl;
    cout << "预处理采样点数：" << abs((dl_first_start.i - dl_first_end.i) / 0x40000000000) << endl;
    double input_x = 0.0, ULP = 0.0, relative = 0.0;
    for (long int i = dl_first_end.i; i <= dl_first_start.i; i += 0x40000000000) {
        DL dl_input;
        dl_input.i = i;
        double input_one = dl_input.d;
        double origin = getDoubleOfOrigin(input_one);
        double ulp = getULP(input_one, origin);
        if (ULP < ulp) {
            ULP = ulp;
            input_x = input_one;
        }
    }
    double origin_relative = getDoubleOfOrigin(input_x);
    relative = getRelativeError(input_x, origin_relative);
    printf("预处理阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
    if (ULP <= 100) {
        cout << "******表达式无显著误差,采用检测方法1检测误差******" << endl;
        double ULP2 = 0.0, input_x2 = 0.0;
        //改变检测策略，不从half开始，增大尾数位
        DL dl_float_start, dl_float_end;
        dl_float_start.d = start;
        dl_float_start.i = dl_float_start.i & 0xFFFFFFFFE0000000;
        dl_float_end.d = end;
        dl_float_end.i = dl_float_end.i & 0xFFFFFFFFE0000000;
        //采样
        /* long int sample_number = 100000;
        long int step = abs((dl_float_end.i - dl_float_start.i)) / sample_number;
        for (long int i = 0; i <= sample_number; ++i) {
            DL dl_input;
            dl_input.i = dl_float_start.i - i * step;
            double origin = getDoubleOfOrigin(dl_input.d);
            double ulp = getULP(dl_input.d, origin);
            if (ULP2 < ulp) {
                ULP2 = ulp;
                input_x2 = dl_input.d;
            }
        } */
        //模拟退火
        vector<double> vec1 = simulated_annealing(dl_float_start.d, dl_float_end.d);
        input_x2 = vec1[0];
        ULP2 = vec1[1];

        printf("第一阶段x = %.8lf, maximumULP = %.2lf\n", input_x2, ULP2);
        DL dl_double_start, dl_double_end;
        dl_double_start.d = input_x2;
        dl_double_end.d = input_x2;
        dl_double_start.i = dl_double_start.i & 0xFFFFFFFFE0000000;
        dl_double_end.i = dl_double_end.i | 0x000000001FFFFFFF;
        printf("%.16lf, %.16lf\n", dl_double_start.d, dl_double_end.d);
        //采样
        /* step = abs((dl_double_end.i - dl_double_start.i)) / sample_number;
        for (long int i = 0; i <= sample_number; ++i) {
            DL dl_input;
            dl_input.i = dl_double_start.i + i * step;
            double origin = getDoubleOfOrigin(dl_input.d);
            double ulp = getULP(dl_input.d, origin);
            if (ULP2 < ulp) {
                ULP2 = ulp;
                input_x2 = dl_input.d;
            }
        } */
        //模拟退火
        vector<double> vec2 = simulated_annealing(dl_double_start.d, dl_double_end.d);
        input_x2 = vec2[0];
        ULP2 = vec2[1];
        origin_relative = getDoubleOfOrigin(input_x2);
        relative = getRelativeError(input_x2, origin_relative);
        printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x2, ULP2, relative, log2(ULP2 + 1));
    } else {
        //半精度遍历完之后得到一个触发最大误差的输入，对该输入周围的数操作更多的尾数位
        cout << "******表达式存在显著误差,采用检测方法2检测误差******" << endl;
        printf("第一阶段x = %.6lf, maximumULP = %.2lf\n", input_x, ULP);
        double input_second = input_x;
        DL dl_float_start, dl_float_end;
        dl_float_start.d = input_second;
        dl_float_end.i = dl_float_start.i | 0x000003FFE0000000;
        cout << "第二阶段采样点数：" << abs((dl_float_end.i - dl_float_start.i)) / (0x20000000) << endl;
        for (long int i = dl_float_end.i; i <= dl_float_start.i; i += 0x20000000) {
            DL dl_input;
            dl_input.i = i;
            double input_two = dl_input.d;
            double origin = getDoubleOfOrigin(input_two);
            double ulp = getULP(input_two, origin);
            if (ULP < ulp) {
                ULP = ulp;
                input_x = input_two;
            }
        }
        origin_relative = getDoubleOfOrigin(input_x);
        relative = getRelativeError(input_x, origin_relative);
        printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
        //此时尾数位增大到float的23位之后，此时获得的出发最大误差的输入为input_x，此时将剩下的全部尾数位不能穷举了，只能采样检测
        //采样数设置为100000个点
        DL dl_double_start, dl_double_end;
        dl_double_start.d = input_x;
        // dl_double_end.i = dl_double_start.i | 0B0000000000000000000000000000000000011111111111111111111111111111;
        dl_double_end.i = dl_double_start.i | 0x000000001FFFFFFF;
        //采样
        /* int sample_number = 100000;
        double step = abs((dl_double_start.d - dl_double_end.d)) / (double)sample_number;
        for (long int i = 0; i <= sample_number; i++) {
            double input_three = dl_double_start.d + i * step;
            double origin = getDoubleOfOrigin(input_three);
            double ulp = getULP(input_three, origin);
            if (ULP < ulp) {
                ULP = ulp;
                input_x = input_three;
            }
        }
        cout << "第三阶段采样点数:100000" << endl; */
        //模拟退火
        vector<double> vec3 = simulated_annealing(dl_double_start.d, dl_double_end.d);
        input_x = vec3[0];
        ULP = vec3[1];
        origin_relative = getDoubleOfOrigin(input_x);
        relative = getRelativeError(input_x, origin_relative);
        printf("第三阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x, ULP, relative, log2(ULP + 1));
    }
}
//start和end跨越0的情况
void DoubleFunction::detectMethod2_2(const double &start, const double &end) {
    DL dl_first_start, dl_first_end;
    dl_first_start.d = start, dl_first_end.d = end;
    dl_first_start.i = dl_first_start.i & 0xFFFFFC0000000000;
    dl_first_end.i = dl_first_end.i & 0x7FFFFC0000000000;
    cout << "检测区间[" << dl_first_start.d << ", " << dl_first_end.d << "]" << endl;
    cout << "预处理采样点数：" << abs((dl_first_end.i - dl_first_start.i) / 0x40000000000) << endl;
    double input_x = 0.0, ULP = 0.0, relative = 0.0, origin_relative = 0.0;
    double ULP2 = 0.0, input_x2 = 0.0;
    for (long int i = dl_first_start.i; i <= dl_first_end.i; i -= 0x40000000000) {
        DL dl_input;
        dl_input.i = i;
        double input_one = dl_input.d;
        double origin = getDoubleOfOrigin(input_one);
        double ulp = getULP(input_one, origin);
        if (ULP < ulp) {
            ULP = ulp;
            input_x = input_one;
        }
    }
    for (long int i = 0; i <= dl_first_end.i; i += 0x40000000000) {
        DL dl_input;
        dl_input.i = i;
        double input_one = dl_input.d;
        double origin = getDoubleOfOrigin(input_one);
        double ulp = getULP(input_one, origin);
        if (ULP < ulp) {
            ULP = ulp;
            input_x = input_one;
        }
    }
    origin_relative = getDoubleOfOrigin(input_x);
    relative = getRelativeError(input_x, origin_relative);
    printf("预处理阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
    if (ULP <= 100) {
        cout << "******表达式无显著误差,采用检测方法1检测误差******" << endl;

        //改变检测策略，不从half开始，增大尾数位
        DL dl_float_start, dl_float_end;
        dl_float_start.d = start;
        dl_float_start.i = dl_float_start.i & 0xFFFFFFFFE0000000;
        dl_float_end.d = end;
        dl_float_end.i = dl_float_end.i & 0x7FFFFFFFE0000000;
        //采样
        /* long int sample_number = 100000;
        long int step = abs((dl_float_end.i - dl_float_start.i) / sample_number);
        for (long int i = 0; i <= sample_number / 2; ++i) {
            DL dl_input;
            dl_input.i = dl_float_start.i - i * step;
            double input_one = dl_input.d;
            double origin1 = getDoubleOfOrigin(input_one);
            double ulp1 = getULP(input_one, origin1);
            double origin2 = getDoubleOfOrigin(-input_one);
            double ulp2 = getULP(-input_one, origin2);
            double ulp = max(ulp1, ulp2);
            double input_test;
            if (ulp == ulp1) input_test = input_one;
            else input_test = -input_one;
            if (ULP2 < ulp) {
                ULP2 = ulp;
                input_x2 = input_test;
            }
        } */
        //模拟退火
        vector<double> vec1 = simulated_annealing(dl_float_start.d, dl_float_end.d);
        input_x2 = vec1[0];
        ULP2 = vec1[1];
        origin_relative = getDoubleOfOrigin(input_x2);
        relative = getRelativeError(input_x2, origin_relative);
        printf("第一阶段x = %.8lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x2, ULP2, relative);
        if (input_x2 <= 0) {
            DL dl_double_start, dl_double_end;
            dl_double_start.d = input_x2;
            dl_double_end.d = input_x2;
            dl_double_start.i = dl_double_start.i & 0xFFFFFFFFE0000000;
            dl_double_end.i = dl_double_end.i | 0x000000001FFFFFFF;
            //采样
            /* step = abs((dl_double_end.i - dl_double_start.i)) / sample_number;
            for (long int i = 0; i <= sample_number; ++i) {
                DL dl_input;
                dl_input.i = dl_double_start.i + i * step;
                double origin = getDoubleOfOrigin(dl_input.d);
                double ulp = getULP(dl_input.d, origin);
                if (ULP2 < ulp) {
                    ULP2 = ulp;
                    input_x2 = dl_input.d;
                }
            } */
            //模拟退火
            vector<double> vec2 = simulated_annealing(dl_double_start.d, dl_double_end.d);
            input_x2 = vec2[0];
            ULP2 = vec2[1];
            origin_relative = getDoubleOfOrigin(input_x2);
            relative = getRelativeError(input_x2, origin_relative);
            printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x2, ULP2, relative, log2(ULP2 + 1));
        } else {
            DL dl_double_start, dl_double_end;
            dl_double_start.d = input_x2;
            dl_double_end.d = input_x2;
            dl_double_start.i = dl_double_start.i & 0x7FFFFFFFE0000000;
            dl_double_end.i = dl_double_end.i | 0x000000001FFFFFFF;
            //采样
            /* step = (dl_double_end.i - dl_double_start.i) / sample_number;
            for (long int i = 0; i <= sample_number; ++i) {
                DL dl_input;
                dl_input.i = dl_double_start.i + i * step;
                double origin = getDoubleOfOrigin(dl_input.d);
                double ulp = getULP(dl_input.d, origin);
                if (ULP2 < ulp) {
                    ULP2 = ulp;
                    input_x2 = dl_input.d;
                }
            } */
            //模拟退火
            vector<double> vec3 = simulated_annealing(dl_double_start.d, dl_double_end.d);
            input_x2 = vec3[0];
            ULP2 = vec3[1];
            origin_relative = getDoubleOfOrigin(input_x2);
            relative = getRelativeError(input_x2, origin_relative);
            printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x2, ULP2, relative, log2(ULP2 + 1));
        }
    } else {
        if (input_x2 <= 0) {
            cout << "******表达式存在显著误差,采用检测方法2检测误差******" << endl;
            printf("第一阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
            double input_second = input_x;
            DL dl_float_start, dl_float_end;
            dl_float_start.d = input_second;
            dl_float_end.i = dl_float_start.i | 0x000003FFE0000000;
            cout << "第二阶段采样点数：" << abs((dl_float_end.i - dl_float_start.i)) / (0x20000000) << endl;
            for (long int i = dl_float_end.i; i <= dl_float_start.i; i += 0x20000000) {
                DL dl_input;
                dl_input.i = i;
                double input_two = dl_input.d;
                double origin = getDoubleOfOrigin(input_two);
                double ulp = getULP(input_two, origin);
                if (ULP < ulp) {
                    ULP = ulp;
                    input_x = input_two;
                }
            }
            origin_relative = getDoubleOfOrigin(input_x);
            relative = getRelativeError(input_x, origin_relative);
            printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
            //此时尾数位增大到float的23位之后，此时获得的出发最大误差的输入为input_x，此时将剩下的全部尾数位不能穷举了，只能采样检测
            //采样数设置为100000个点
            DL dl_double_start, dl_double_end;
            dl_double_start.d = input_x;
            // dl_double_end.i = dl_double_start.i | 0B0000000000000000000000000000000000011111111111111111111111111111;
            dl_double_end.i = dl_double_start.i | 0x000000001FFFFFFF;
            //采样
            /* int sample_number = 100000;
            double step = abs((dl_double_start.d - dl_double_end.d)) / (double)sample_number;
            for (long int i = 0; i <= sample_number; i++) {
                double input_three = dl_double_start.d + i * step;
                double origin = getDoubleOfOrigin(input_three);
                double ulp = getULP(input_three, origin);
                if (ULP < ulp) {
                    ULP = ulp;
                    input_x = input_three;
                }
            }
            cout << "第三阶段采样点数:100000" << endl; */
            //模拟退火
            vector<double> vec4 = simulated_annealing(dl_double_start.d, dl_double_end.d);
            input_x = vec4[0];
            ULP = vec4[1];
            origin_relative = getDoubleOfOrigin(input_x);
            relative = getRelativeError(input_x, origin_relative);
            printf("第三阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x, ULP, relative, log2(ULP + 1));
        } else {    
            cout << "******表达式存在显著误差,采用检测方法2检测误差******" << endl;
            printf("第一阶段x = %.6lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
            double input_second = input_x;
            DL dl_float_start, dl_float_end;
            dl_float_start.d = input_second;
            // dl_float_end.i = dl_float_start.i | 0B0000000000000000000000111111111111100000000000000000000000000000;
            dl_float_end.i = dl_float_start.i | 0x000003FFE0000000;
            cout << "第二阶段采样点数：" << (dl_float_end.i - dl_float_start.i) / (0x20000000) << endl;
            for (long int i = dl_float_start.i; i <= dl_float_end.i; i += 0x20000000) {
                DL dl_input;
                dl_input.i = i;
                double input_two = dl_input.d;
                double origin = getDoubleOfOrigin(input_two);
                double ulp = getULP(input_two, origin);
                if (ULP < ulp) {
                    ULP = ulp;
                    input_x = input_two;
                }
            }
            origin_relative = getDoubleOfOrigin(input_x);
            relative = getRelativeError(input_x, origin_relative);
            printf("第二阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e\n", input_x, ULP, relative);
            //此时尾数位增大到float的23位之后，此时获得的出发最大误差的输入为input_x，此时将剩下的全部尾数位不能穷举了，只能采样检测
            //采样数设置为100000个点
            DL dl_double_start, dl_double_end;
            dl_double_start.d = input_x;
            // dl_double_end.i = dl_double_start.i | 0B0000000000000000000000000000000000011111111111111111111111111111;
            dl_double_end.i = dl_double_start.i | 0x000000001FFFFFFF;
            //采样
            /* int sample_number = 100000;
            double step = (dl_double_start.d - dl_double_end.d) / (double)sample_number;
            for (long int i = 0; i <= sample_number; i++) {
                double input_three = dl_double_start.d + i * step;
                double origin = getDoubleOfOrigin(input_three);
                double ulp = getULP(input_three, origin);
                if (ULP < ulp) {
                    ULP = ulp;
                    input_x = input_three;
                }
            }
            cout << "第三阶段采样点数:100000" << endl; */
            //模拟退火
            vector<double> vec5 = simulated_annealing(dl_double_start.d, dl_double_end.d);
            input_x = vec5[0];
            ULP = vec5[1];
            origin_relative = getDoubleOfOrigin(input_x);
            relative = getRelativeError(input_x, origin_relative);
            printf("第三阶段x = %.16lf, maximumULP = %.2lf, maximumRelative = %e, BitsError = %.1lf\n", input_x, ULP, relative, log2(ULP + 1));
        }   
    }
}
//模拟退火算法
vector<double> DoubleFunction::simulated_annealing(const double &start, const double &end) {
    vector<double> result(2);
    int num = 2000;//迭代次数10000->2000
    double k = 0.1;
    double r = 0.9;//控制降温快慢
    double T = 3000;//系统的温度，系统初始化温度设置为一个相对较高的温度
    double T_min = 20;//温度的下限，若温度T达到T_min，则停止搜索
    //生成[start,end]范围内的随机浮点数
    // random_device rd;
    // default_random_engine eng(rd());
    // uniform_real_distribution<> distr(start, end);
    // double dbNumber = distr(eng);
    //伪随机数
    double dbNumber = start + static_cast<double> (rand()) / (static_cast<double> (RAND_MAX / (end - start)));

    double origin = getDoubleOfOrigin(dbNumber);
    double ULP = getULP(dbNumber, origin);//目前函数值
    double maximum = ULP;
    double inputx = 0.0;
    double dE, current;
    while (T > T_min) {
        for (int i = 0; i < num; ++i) {
            // random_device rd;
            // default_random_engine eng(rd());
            // uniform_real_distribution<> distr(start, end);
            // double dbNumber = distr(eng);
            double dbNumber2 = start + static_cast<double> (rand()) / (static_cast<double> (RAND_MAX / (end - start)));
            double origin2 = getDoubleOfOrigin(dbNumber2);
            current = getULP(dbNumber2, origin2);
            // cout << "dbNumber = " << setprecision(17) << dbNumber2 << " ulp = " << current << endl;
            // cout << "current = " << setprecision(16) << current << endl;
            dE = ULP - current;
            double dEtmp = maximum - current;
            if (dEtmp < 0) {
                maximum = current;
                inputx = dbNumber2;
            }
            if (dE < 0) {
                ULP = current;
            } else {
                double factor = static_cast<double> (rand()) / (static_cast<double> (RAND_MAX / 1.0));
                if (exp(-dE / (T * k)) > factor) ULP = current;
            }
        }
        T = r * T;//降温退火，0<r<1。r越大，降温越慢；r越小，降温越快
    }
    result[0] = inputx;
    result[1] = maximum;
    return result;
}

